Atelier
Sana & Hana
2023-03-29
Importation de la base:
library(readxl)
Book1 <- read_excel("C:/Users/sanat/OneDrive/Bureau/Hello/Atelier/Book1.xlsx")
BD <- Book1Chargement des librairies:
library(fpp3)## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.2.0 ✔ tsibble 1.1.3
## ✔ dplyr 1.1.0 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.0 ✔ feasts 0.3.1
## ✔ lubridate 1.9.2 ✔ fable 0.3.3
## ✔ ggplot2 3.4.2 ✔ fabletools 0.3.2## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()library(zoo)##
## Attachement du package : 'zoo'## L'objet suivant est masqué depuis 'package:tsibble':
##
## index## Les objets suivants sont masqués depuis 'package:base':
##
## as.Date, as.Date.numericlibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attachement du package : 'forecast'## Les objets suivants sont masqués depuis 'package:fabletools':
##
## accuracy, forecastlibrary(tseries)
library(astsa)##
## Attachement du package : 'astsa'## L'objet suivant est masqué depuis 'package:forecast':
##
## gaslibrary(dplyr)
library(lubridate)
library(tsibble)
library(tsibble)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ purrr 1.0.1 ✔ stringr 1.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(ggplot2)
library(gridExtra)##
## Attachement du package : 'gridExtra'
##
## L'objet suivant est masqué depuis 'package:dplyr':
##
## combinelibrary(grid)
library(ggthemes)
library(corrplot)## corrplot 0.92 loadedlibrary(lmtest)
library(sandwich)
library(TSA)## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attachement du package : 'TSA'
##
## L'objet suivant est masqué depuis 'package:readr':
##
## spec
##
## Les objets suivants sont masqués depuis 'package:stats':
##
## acf, arima
##
## L'objet suivant est masqué depuis 'package:utils':
##
## tarlibrary(car)## Le chargement a nécessité le package : carData
##
## Attachement du package : 'car'
##
## L'objet suivant est masqué depuis 'package:purrr':
##
## some
##
## L'objet suivant est masqué depuis 'package:dplyr':
##
## recodelibrary(aTSA)##
## Attachement du package : 'aTSA'
##
## Les objets suivants sont masqués depuis 'package:tseries':
##
## adf.test, kpss.test, pp.test
##
## L'objet suivant est masqué depuis 'package:forecast':
##
## forecast
##
## Les objets suivants sont masqués depuis 'package:fabletools':
##
## estimate, forecast
##
## L'objet suivant est masqué depuis 'package:graphics':
##
## identifylibrary(forecast)
library(astsa)
library(urca)
library(seastests)
library(TSstudio)
library("cowplot")##
## Attachement du package : 'cowplot'
##
## L'objet suivant est masqué depuis 'package:TSstudio':
##
## plot_grid
##
## L'objet suivant est masqué depuis 'package:ggthemes':
##
## theme_map
##
## L'objet suivant est masqué depuis 'package:lubridate':
##
## stamplibrary(stats)
library(uroot)Transformer le vecteur BD$Nombre en un objet de type ts.
La commande ts permet de convertir un vecteur en une série temporelle.
ts1 <- ts(BD$Nombre,start=c(2007,1), end=c(2016,6), frequency = 12)
ts1## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2007 353761 325955 452936 717818 1435020 1966038 2744558 3229190 2468861
## 2008 390173 383587 464958 689467 1461643 2089183 2617998 3171273 2345464
## 2009 323236 314462 394722 680024 1309381 1892282 2512378 3003760 2266738
## 2010 274952 296528 408930 563415 1247893 1863445 2783843 3155878 2298412
## 2011 362517 333197 412673 669391 1342966 2183411 2984058 3352934 2588707
## 2012 352320 266066 360171 625837 1179714 2043692 2863095 3269202 2522887
## 2013 365604 295651 362098 557496 1467134 2372846 3263920 3885717 2963438
## 2014 404292 325952 456654 727963 1651703 2697469 4222873 4856356 3643695
## 2015 606140 509188 613091 934237 1870169 3032869 4408555 4993464 3649700
## 2016 559060 433188 628721 901585 1952667 2970996
## Oct Nov Dec
## 2007 1473607 517954 468235
## 2008 1423374 487813 402002
## 2009 1344580 450735 410632
## 2010 1316460 458142 339591
## 2011 1428960 423107 345325
## 2012 1310986 392800 330846
## 2013 1547626 435303 402743
## 2014 1840599 668173 537728
## 2015 1852679 641458 487900
## 2016Données représentant le nombre mensuel des visiteurs vers la Grèce de janvier 2007 à juin 2016. Chaque ligne de données représente une année soit douze valeurs Lors de la conversion du vecteur en série temporelle, on a start = c(2007,1) ,end=c(2016,6) , frequency=12 (cad mensuelle)
str(ts1)## Time-Series [1:114] from 2007 to 2016: 353761 325955 452936 717818 1435020 ...Cette sortie nous indique que “ts1” est un objet de série chronologique avec 114 observations (mois).
summary(ts1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 266066 425627 917911 1446797 2333701 4993464La fonction de synthèse a été utilisée pour obtenir la valeur minimale, le premier quartile, la médiane, la moyenne, le troisième quartile et la valeur maximale de la série. La série va d’une valeur minimale de 1 à une valeur maximale de 4 993 464, avec une valeur médiane de 917 911. Sans contexte supplémentaire, il est difficile de fournir une interprétation plus précise de ces données.
Représenter graphiquement la série.
La représentation de la série :Ici on va diviser la fenêtre en deux parties pour bien comprendre le concept de la sérietemporelle , donc dans la première partie on va représenter le vecteur BD$Nombre , et la deuxièmepartie on va représenter la série temporelle ts1:
par (mfrow=c(1,2))
plot(BD$Nombre)
plot.ts(ts1)
==> La deuxième partie est la représentation de la série temporelle
(les données en fonction du temps). La saisonnalité devient très
apparente dans cette nouvelle série , il y a un comportement saisonnier
mais pas de tendance bien évidemment puisque la variations de cette
série a une amplitude qui semble augmenter au cours du temps.
Variance d’une série temporelle:
variance <- rollapply(ts1, width = 12, FUN = var, align = "right")
print(variance)## Jan Feb Mar Apr May
## 2007
## 2008 1.070375e+12 1.059930e+12 1.057972e+12 1.061324e+12 1.061781e+12
## 2009 1.029116e+12 1.041304e+12 1.052582e+12 1.053655e+12 1.051369e+12
## 2010 9.410865e+11 9.441244e+11 9.419669e+11 9.549211e+11 9.543248e+11
## 2011 1.072328e+12 1.066031e+12 1.065452e+12 1.052941e+12 1.053310e+12
## 2012 1.286617e+12 1.299624e+12 1.308920e+12 1.314532e+12 1.317095e+12
## 2013 1.220537e+12 1.215079e+12 1.214751e+12 1.223479e+12 1.224538e+12
## 2014 1.716499e+12 1.709960e+12 1.691159e+12 1.664154e+12 1.665180e+12
## 2015 2.621147e+12 2.573072e+12 2.534963e+12 2.495255e+12 2.489432e+12
## 2016 2.732948e+12 2.753515e+12 2.749718e+12 2.755882e+12 2.755178e+12
## Jun Jul Aug Sep Oct
## 2007
## 2008 1.076729e+12 1.046322e+12 1.026861e+12 1.003022e+12 1.002007e+12
## 2009 1.026224e+12 1.001462e+12 9.459356e+11 9.308778e+11 9.289236e+11
## 2010 9.508836e+11 1.020797e+12 1.071425e+12 1.077331e+12 1.076956e+12
## 2011 1.095898e+12 1.153083e+12 1.222045e+12 1.279769e+12 1.279892e+12
## 2012 1.297319e+12 1.262143e+12 1.231766e+12 1.216850e+12 1.215387e+12
## 2013 1.277169e+12 1.401370e+12 1.645265e+12 1.749193e+12 1.747274e+12
## 2014 1.723324e+12 2.096422e+12 2.570644e+12 2.762459e+12 2.757178e+12
## 2015 2.546422e+12 2.626211e+12 2.699979e+12 2.701807e+12 2.701531e+12
## 2016 2.743448e+12
## Nov Dec
## 2007 1.076834e+12
## 2008 1.006562e+12 1.017338e+12
## 2009 9.341378e+11 9.328273e+11
## 2010 1.075876e+12 1.087223e+12
## 2011 1.285812e+12 1.284742e+12
## 2012 1.220279e+12 1.222794e+12
## 2013 1.738994e+12 1.724307e+12
## 2014 2.703684e+12 2.670299e+12
## 2015 2.707928e+12 2.721118e+12
## 2016plot(variance)
Notre série montre une variance croissante au fil des années, cela peut
indiquer une instabilité ou une non-stationnarité dans les données. Cela
signifie que les propriétés statistiques de la série, telles que la
moyenne et la variance, changent avec le temps.
ACF:
autocorr = acf(ts1, lag.max=25)
La fonction d’autocorrélation est périodique,ce qui indique une
périodicité dans la série temporelle. La ligne pointillée bleue indique
le niveau en-dessous duquel la corrélation n’est plus statistiquement
significative. Les valeurs ≥ 0.7 sont fortement correles Les valeurs <
0.7 sont faiblement correles Les valeurs a l’interieur de l’interval de
confiance ne sont pas correles. On trouve que chaque début d’année les
données ont une correlation positive forte puis diminue a cause de la
perte de données Les traits pointillés horizentaux représentent
l’intervalle de confiance ; les valeurs à l’intérieurde cet intervalle
indique que la corrélation et nulle donc les données sont indépendants ,
nesont pas corrélées (liées)
Vérifier si la série est corrélée:
La vérification de la corrélation entre les observations d’une série chronologique est importante car cela peut fournir des informations sur la structure de la série temporelle. Si les observations sont corrélées, cela indique qu’il existe une certaine dépendance entre les observations et que les observations passées peuvent être utilisées pour prédire les observations futures. Cela est particulièrement important pour la modélisation et la prévision de la série chronologique, car cela permet l’utilisation de modèles qui prennent en compte la dépendance temporelle des données. Pour vérifier la corrélation de notre série, on a utilisé le Box.test qui effectue un test statistique de la boîte de Ljung-Box ou si le p-value du test est inférieur à un certain seuil (généralement 0,05), on rejette l’hypothèse nulle et on conclut que la série chronologique est autocorrélée jusqu’à un certain lag. Plus précisément, le test de la boîte de Ljung-Box est utilisé pour tester l’hypothèse nulle que les corrélations croisées d’ordre k dans la série chronologique sont nulles pour tous les k jusqu’à un certain nombre de retards spécifié. Dans ce cas, le test est effectué sur les autocorrélations de la série, qui sont stockées dans la variable “autocorr$acf”.
Box.test(autocorr$acf, lag = 12, type = "Ljung-Box")##
## Box-Ljung test
##
## data: autocorr$acf
## X-squared = 115.44, df = 12, p-value < 2.2e-16le test de Box-Ljung a donné une statistique de test X² de 115.44 pour 12 degrés de liberté et une p-value de 2,2e-16 < 0.05 : seuil de significativité. On peut rejeter l’hypothèse nulle donc on peut conclure que les autocorrélations jusqu’à un certain lag ne sont pas nulles et que la série chronologique est autocorrélée.
Décomposition:
Commencer avec la fonction decompose(): Décomposer la série en faisant varier à chaque fois le type de la décomposition avec les options type=“add” ou “mult”, qui permettent de spécifier si on souhaite utiliser un modèle additif (par défaut) ou multiplicatif. Decompose() :Décomposer une série chronologique en composantes saisonnières, tendancielles et irrégulières en utilisant des moyennes mobiles. Traite avec un composant saisonnier additif ou multiplicatif. Le modèle additif utilisé est: Y [t] = T [t] + S [t] + e [t] Le modèle multiplicatif utilisé est: Y [t] = T [t] * S [t] * e [t] decompose() utilise la technique de la moyenne mobile pour l’estimation de tendance, d’où la présence de NA au niveau de TREND ..
Le modèle additif :
plot(decompose(ts1,"additive"))
Le modèle multiplicatif :
plot(decompose(ts1,"multiplicative"))
Analyser les résidus de chaque type de décomposition
residusadd= (decompose(ts1, type="additive"))$random
residusmult= (decompose(ts1, type="mult"))$random
table_residus <- data.frame(residusadd, residusmult)
names(table_residus) <- c("Additive", "Multiplicative")
print(table_residus)## Additive Multiplicative
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## 7 -299593.786 0.9448772
## 8 -319266.268 0.9559690
## 9 -173423.694 0.9696584
## 10 79115.163 1.0500113
## 11 133824.339 1.1227853
## 12 168301.102 1.2131464
## 13 72921.096 1.0724121
## 14 136202.591 1.2014925
## 15 146884.221 1.1564572
## 16 139935.716 1.1019677
## 17 157731.591 1.1107843
## 18 -40541.404 1.0149835
## 19 -402931.078 0.9171093
## 20 -344372.810 0.9621805
## 21 -255301.277 0.9503238
## 22 73041.163 1.0484153
## 23 154508.047 1.0986969
## 24 173681.602 1.0994375
## 25 90060.346 0.9469819
## 26 154847.841 1.0533496
## 27 167123.679 1.0525907
## 28 220296.924 1.1651122
## 29 96752.883 1.0681633
## 30 -148989.362 0.9848197
## 31 -423994.536 0.9401298
## 32 -430239.435 0.9715022
## 33 -258032.360 0.9748501
## 34 71188.829 1.0521937
## 35 195054.714 1.0795002
## 36 249151.352 1.1844901
## 37 85901.554 0.8343836
## 38 152009.216 1.0052782
## 39 178509.096 1.0880299
## 40 94153.924 0.9580003
## 41 21765.924 1.0070638
## 42 -189859.237 0.9605028
## 43 -166903.161 1.0297741
## 44 -300430.560 1.0026214
## 45 -250506.569 0.9695448
## 46 10082.288 1.0033400
## 47 153677.422 1.0552268
## 48 108269.227 0.9265907
## 49 92060.721 1.0319086
## 50 108368.716 1.0613634
## 51 89294.304 1.0223788
## 52 90537.091 1.0469259
## 53 3155.341 0.9938499
## 54 14992.638 1.0307736
## 55 -80927.828 1.0116906
## 56 -209215.810 0.9826466
## 57 -59384.277 1.0123886
## 58 31983.538 1.0163330
## 59 45037.630 0.9213069
## 60 70315.519 0.9114031
## 61 70711.971 0.9944370
## 62 55167.883 0.8564982
## 63 77260.096 0.9198839
## 64 111892.091 1.0276125
## 65 -85781.576 0.9231136
## 66 -49766.071 1.0207154
## 67 -127066.703 1.0268223
## 68 -223131.893 1.0097417
## 69 -61686.069 1.0349750
## 70 -23707.337 0.9773043
## 71 59268.547 0.8845462
## 72 62060.061 0.8772788
## 73 48941.971 1.0053189
## 74 -1219.242 0.8933596
## 75 -57060.779 0.8376478
## 76 -128571.159 0.8084979
## 77 11706.841 1.0012382
## 78 82823.596 1.0287994
## 79 72536.464 1.0130097
## 80 191072.940 1.0380522
## 81 172665.390 1.0489178
## 82 -7076.754 0.9860625
## 83 -123902.745 0.8355596
## 84 -87242.981 0.9155234
## 85 -156636.612 0.9423543
## 86 -253195.326 0.8198164
## 87 -269524.779 0.8713400
## 88 -277459.117 0.8677062
## 89 -133358.284 0.9224980
## 90 67251.846 0.9521886
## 91 681867.714 1.0622116
## 92 798919.565 1.0443841
## 93 481179.348 1.0339609
## 94 -89917.129 0.9395205
## 95 -269750.453 1.0279171
## 96 -332837.106 0.9814097
## 97 -303597.154 1.1604155
## 98 -351817.784 1.0970541
## 99 -332121.946 1.0398837
## 100 -250421.576 1.0123891
## 101 -71608.826 0.9615009
## 102 264451.888 0.9954290
## 103 747422.297 1.0411138
## 104 837073.649 1.0196407
## 105 404898.890 0.9921190
## 106 -144300.379 0.9135580
## 107 -347308.120 0.9611994
## 108 -411289.398 0.8774588
## 109 NA NA
## 110 NA NA
## 111 NA NA
## 112 NA NA
## 113 NA NA
## 114 NA NApar (mfrow=c(1,2))
hist(residusadd)
hist(residusmult)
==>On remarque que les residusadd ne suivent pas la loi normale
==>Le modèle multiplicatif est meilleur quele modèle additif car les
résidus du modèle multiplicatif présentent moins de fluctuations que le
modèle additif. De meme l’histogramme des résidus nous conduit a dire
que les données du modèle multiplicatif sont plus homogènes que celles
du modèle additif.
En outre, la normalité des résidus est également importante pour l’interprétation des résultats. Si les résidus ne sont pas normalement distribués, il peut être difficile d’interpréter correctement les coefficients de régression ou les intervalles de confiance associés. En revanche, si les résidus suivent une distribution normale, cela permet une interprétation plus fiable des résultats du modèle.
Le test de Shapiro-Wilk est un test statistique qui permet de vérifier si un échantillon de données suit ou non une distribution normale. Sous R, la fonction shapiro.test() permet de réaliser ce test. Le résultat du test inclura la statistique de test W, la valeur de p et une indication sur la signification du résultat. Si la valeur de p est inférieure à 0,05, alors on peut rejeter l’hypothèse nulle selon laquelle les données suivent une distribution normale. Si la valeur de p est supérieure à 0,05, on ne peut pas rejeter l’hypothèse nulle et on peut considérer que les données suivent une distribution normale.
shapiro.test(residusmult)##
## Shapiro-Wilk normality test
##
## data: residusmult
## W = 0.98609, p-value = 0.3656Le résultat du test de Shapiro-Wilk indique que la statistique de test W est égale à 0,98609 et que la valeur de p est égale à 0,3656. La valeur de p étant supérieure à 0,05, on ne peut pas rejeter l’hypothèse nulle selon laquelle les résidus multiplicatifs suivent une distribution normale. Autrement dit, il n’y a pas suffisamment de preuves pour dire que les données ne suivent pas une distribution normale. Cependant, il est important de garder à l’esprit que les tests de normalité sont souvent sensibles à la taille de l’échantillon. Si l’échantillon est petit, le test peut manquer de puissance pour détecter des écarts par rapport à une distribution normale. Dans ce cas, il peut être utile d’utiliser d’autres méthodes de diagnostic pour évaluer la normalité des données.
Observation de la saisonnalité:
On va visualiser les composantes saisonnières de la série temporelle ts1 à l’aide d’un graphique en boîte (ou boxplot en anglais). La fonction “ts_seasonal” calcule les moyennes des observations de la série pour chaque période saisonnière et trace les boîtes qui montrent la variabilité des observations dans chaque période saisonnière. La visualisation permet ainsi d’observer les variations saisonnières de la série. Ensuite, On va calculer la matrice de corrélation de la série ts1. Et comme ts1 est une série temporelle univariée, la fonction renverra simplement la corrélation entre chaque observation et son observation précédente.
ts_seasonal(ts1,type="box")ts_cor(ts1)Graphiquement, esa boites à moustache prouve l’existence de la saisonnalité. S’il n’y avait pas d’effet saisonnier, ces 12 boites se ressembleraient ! De plus, la présence de seasonal lag 12 dans l’ACF indique la présence de la saisonnalité.
Teste de saisonnalité:
a fonction isSeasonal est une fonction R qui permet de déterminer si une série temporelle a une composante saisonnière significative. La fonction utilise une approche basée sur les filtres de Hodrick-Prescott pour décomposer la série temporelle en une tendance et une composante saisonnière. Ensuite, elle calcule une statistique de test basée sur la variance de la composante saisonnière pour déterminer si la série a une composante saisonnière significative: 1- Calcule la variance de la série temporelle originale. 2- Calcule la variance de la série temporelle différée d’une saison en utilisant la fonction diff() avec un décalage de la période de saisonnalité. 3- Compare les variances en calculant le rapport entre la variance de la série temporelle différée d’une saison et la variance de la série temporelle originale. 4- Si le rapport de variance est significativement inférieur à 1, la série est considérée comme saisonnière. Si la statistique de test est supérieure à une valeur seuil, la fonction renvoie TRUE, indiquant que la série temporelle est saisonnière. Sinon, elle renvoie FALSE, indiquant que la série n’a pas de composante saisonnière significative.
isSeasonal(ts1) #renvoie une valeur logique## [1] TRUELe test de Canova et Hansen est un test de racine unitaire saisonnière. Il permet de tester si une série chronologique saisonnière suit un processus ARIMA(p,1,q)(P,1,Q)_s avec une racine unitaire. Le test est basé sur la statistique de test suivante : CH = n*(SSE1 - SSE2)/SSE2 où SSE1 et SSE2 sont les sommes des carrés des erreurs de régression pour les modèles ARIMA(p,1,q)(P,1,Q)_s avec et sans racine unitaire, respectivement, n est la taille de l’échantillon et s est la période saisonnière. La statistique de test suit une distribution chi-carré avec p+P+q+Q degrés de liberté. Sous R, le test de Canova et Hansen peut être effectué à l’aide de la fonction ch.test(). Cette fonction prend en entrée une série chronologique et retourne les résultats du test ainsi que des informations sur le modèle ARIMA estimé avec et sans racine unitaire. Voici un exemple de code :
ch.test(ts1)##
## Canova and Hansen test for seasonal stability
##
## data: ts1
##
## statistic pvalue
## Jan 0.589 0.0127 *
## Feb 0.3277 0.1209
## Mar 0.4101 0.0624 .
## Apr 0.3977 0.0694 .
## May 0.586 0.0131 *
## Jun 0.804 0.0012 **
## Jul 0.6225 0.0089 **
## Aug 0.6079 0.0104 *
## Sep 0.625 0.0087 **
## Oct 0.4744 0.0356 *
## Nov 0.2825 0.1706
## Dec 0.2333 0.2493
## joint 1.7768 0.4102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test type: seasonal dummies
## NW covariance matrix lag order: 12
## First order lag: no
## Other regressors: no
## P-values: based on response surface regressionsLe test de Canova et Hansen pour la stabilité saisonnière évalue si les coefficients de la saisonnalité sont constants au fil du temps ou s’ils varient. La sortie montre les résultats pour chaque mois de l’année, ainsi que la statistique de test jointe pour l’ensemble des mois. Pour chaque mois, on peut voir la statistique de test et la valeur p correspondante. Si la valeur p est inférieure à un niveau de signification spécifié, cela indique que la stabilité saisonnière pour ce mois est rejetée au profit de l’hypothèse alternative de non-stationnarité saisonnière. Dans cet exemple, la stabilité saisonnière est rejetée pour les mois de janvier, mai, juin, juillet, août, septembre et octobre, car leur p-value est inférieure à 0,05. Pour les autres mois, la stabilité saisonnière n’est pas rejetée. La statistique de test jointe pour l’ensemble des mois n’est pas significative, avec une p-value de 0,4102, ce qui suggère qu’il n’y a pas de preuve de non-stationnarité saisonnière globale dans la série chronologique étudiée. Étant donné que la valeur de p est supérieure au niveau de signification de 0,05, on ne peut pas rejeter l’hypothèse nulle selon laquelle toutes les périodes ont une racine unitaire commune. Par conséquent, on peut conclure que la série chronologique a des racines unitaires saisonnières.
Plot de la composante saisonnière:
d <- decompose(ts1,"multiplicative")
d$seasonal## Jan Feb Mar Apr May Jun Jul
## 2007 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2008 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2009 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2010 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2011 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2012 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2013 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2014 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2015 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205 2.1553155
## 2016 0.2675508 0.2361105 0.2990129 0.4678358 0.9863914 1.5476205
## Aug Sep Oct Nov Dec
## 2007 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2008 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2009 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2010 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2011 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2012 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2013 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2014 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2015 2.4992094 1.8797460 1.0366383 0.3407666 0.2838023
## 2016plot(d$seasonal)
Ce plot montre que la saisonnalité varie au fil du temps ce qui nous
empeche pour désaisonnaliser notre série d’utiliser la méthode de la
moyenne mobile centrée qui consiste à soustraire la moyenne mobile
centrée de la série chronologique et qui est simple mais elle ne
fonctionne pas bien si la saisonnalité varie au fil du temps. Et de meme
pour la méthode de la régression saisonnière :qui utilise une régression
linéaire pour estimer les coefficients saisonniers et qui peut être
utilisée que si la saisonnalité varie de manière linéaire. Par contre,
on remarque que la composante saisonnière est stable car les motifs
saisonniers sont cohérents d’une année à l’autre. Dans cet exemple, les
valeurs saisonnières semblent être les mêmes pour chaque mois d’une
année à l’autre, ce qui peut indiquer que la composante saisonnière est
stable. Pour cela, nous allons effectuer la méthode STL, cette méthode
utilise une décomposition en sous-séries saisonnières et tendancielles.
Elle est particulièrement utile pour les séries présentant une
saisonnalité non régulière ou des tendances non linéaires. Appliquons la
méthode STL en utilisant la fonction stl() du package “forecast”. Cette
fonction prend en entrée la série chronologique, le nombre de périodes
dans une saisonnalité et le paramètre s.window qui permet de spécifier
la longueur de la fenêtre de lissage pour la saisonnalité.
Désaisonnalisation:
La méthode de désaisonnalisation basée sur la décomposition STL (Seasonal and Trend decomposition using Loess) est une méthode non-paramétrique qui permet de décomposer une série temporelle en trois composantes: la tendance, la saisonnalité et la composante résiduelle. La méthode commence par déterminer une estimation de la tendance en utilisant une moyenne mobile. Ensuite, la tendance est soustraite de la série originale pour obtenir une série corrigée de tendance. La série corrigée de tendance est ensuite lissée pour extraire la saisonnalité en utilisant une méthode de lissage de Loess. Enfin, la saisonnalité est soustraite de la série corrigée de tendance pour obtenir la composante résiduelle.
stl <- stl(ts1, s.window=12)
stl## Call:
## stl(x = ts1, s.window = 12)
##
## Components
## seasonal trend remainder
## Jan 2007 -988674.0806 1331039 11396.415
## Feb 2007 -1019940.0394 1333430 12464.684
## Mar 2007 -927460.3049 1335822 44574.260
## Apr 2007 -690435.5797 1338214 70039.845
## May 2007 21561.5973 1341653 71805.645
## Jun 2007 727792.8802 1345092 -106846.661
## Jul 2007 1520405.2735 1348531 -124378.078
## Aug 2007 1994107.1909 1351080 -115996.860
## Sep 2007 1167998.8493 1353629 -52766.383
## Oct 2007 67293.9154 1356177 50135.687
## Nov 2007 -894867.3186 1355721 57100.615
## Dec 2007 -972676.9043 1355264 85647.895
## Jan 2008 -993045.6335 1354807 28411.318
## Feb 2008 -1027083.0434 1349758 60912.193
## Mar 2008 -935412.6671 1344708 55662.282
## Apr 2008 -698870.7867 1339659 48678.867
## May 2008 16378.1968 1333702 111562.788
## Jun 2008 734334.6806 1327745 27103.209
## Jul 2008 1539010.1852 1321788 -242800.391
## Aug 2008 2014253.0335 1317026 -160005.848
## Sep 2008 1181023.8691 1312263 -147823.291
## Oct 2008 63170.2375 1307501 52702.732
## Nov 2008 -904542.9903 1297762 94593.588
## Dec 2008 -983459.4568 1288024 97437.683
## Jan 2009 -997574.3419 1278285 42525.197
## Feb 2009 -1034231.7325 1268946 79748.156
## Mar 2009 -943219.2440 1259606 78335.236
## Apr 2009 -706986.6099 1250266 136744.170
## May 2009 11687.7788 1245568 52124.900
## Jun 2009 741502.9300 1240870 -90091.132
## Jul 2009 1558375.0123 1236172 -282169.095
## Aug 2009 2035179.5183 1235570 -266989.638
## Sep 2009 1194850.2578 1234968 -163080.415
## Oct 2009 59755.8406 1234366 50457.965
## Nov 2009 -913601.4691 1236221 128115.069
## Dec 2009 -993809.0396 1238077 166364.434
## Jan 2010 -1007690.9168 1239932 42711.106
## Feb 2010 -1049858.1234 1243691 102695.123
## Mar 2010 -960905.8927 1247450 122385.702
## Apr 2010 -726131.8341 1251209 38337.453
## May 2010 -565.6678 1252130 -3670.915
## Jun 2010 756524.3982 1253050 -146129.182
## Jul 2010 1601383.2020 1253970 -71510.187
## Aug 2010 2080844.3353 1260403 -185369.395
## Sep 2010 1222981.7978 1266836 -191405.931
## Oct 2010 49796.2764 1273269 -6605.484
## Nov 2010 -934948.3424 1286702 106388.606
## Dec 2010 -1018246.5490 1300134 57703.285
## Jan 2011 -1018262.7516 1313567 67212.959
## Feb 2011 -1065439.4348 1326488 72148.546
## Mar 2011 -978047.1225 1339409 51311.138
## Apr 2011 -744129.8776 1352330 61190.796
## May 2011 -11070.1718 1356414 -2378.179
## Jun 2011 773837.6980 1360499 49074.683
## Jul 2011 1647226.1124 1364583 -27751.001
## Aug 2011 2129673.9174 1362577 -139317.262
## Sep 2011 1254608.1472 1360572 -26472.948
## Oct 2011 43393.4986 1358566 27000.244
## Nov 2011 -952676.4525 1348805 26978.912
## Dec 2011 -1039293.1824 1339043 45575.359
## Jan 2012 -1052559.7498 1329281 75598.643
## Feb 2012 -1110797.9999 1318868 57995.555
## Mar 2012 -1018256.3889 1308456 69971.607
## Apr 2012 -779899.1289 1298043 107693.009
## May 2012 -15097.4004 1295754 -100943.016
## Jun 2012 819748.2630 1293466 -69521.976
## Jul 2012 1715954.3855 1291177 -144036.395
## Aug 2012 2207332.7951 1297653 -235784.217
## Sep 2012 1299196.8311 1304130 -80439.665
## Oct 2012 33426.8031 1310606 -33047.050
## Nov 2012 -981203.3877 1331141 42862.850
## Dec 2012 -1072242.3274 1351675 51413.499
## Jan 2013 -1088852.3003 1372209 82247.181
## Feb 2013 -1158269.6734 1399055 54865.375
## Mar 2013 -1060696.3198 1425901 -3107.157
## Apr 2013 -818059.6905 1452748 -77191.964
## May 2013 -21676.5854 1466369 22441.543
## Jun 2013 862963.2530 1479990 29892.316
## Jul 2013 1781843.4651 1493612 -11535.284
## Aug 2013 2282060.3186 1503533 100123.357
## Sep 2013 1340762.0000 1513455 109221.170
## Oct 2013 20447.9266 1523376 3801.738
## Nov 2013 -1012731.1700 1556932 -108897.540
## Dec 2013 -1107956.5380 1590487 -79787.547
## Jan 2014 -1103717.8630 1624042 -116032.597
## Feb 2014 -1176818.0535 1671475 -168704.929
## Mar 2014 -1076576.0568 1718908 -185677.449
## Apr 2014 -831646.8927 1766340 -206730.135
## May 2014 -22162.5602 1794373 -120507.288
## Jun 2014 878688.7069 1822406 -3625.374
## Jul 2014 1806733.7213 1850438 565700.792
## Aug 2014 2310707.3362 1865702 679946.436
## Sep 2014 1358307.2550 1880966 404421.777
## Oct 2014 16493.0220 1896230 -72123.731
## Nov 2014 -1024219.0788 1911556 -219164.417
## Dec 2014 -1120927.1966 1926883 -268228.085
## Jan 2015 -1116738.5801 1942210 -219331.489
## Feb 2015 -1193880.2131 1949479 -246410.422
## Mar 2015 -1091328.1985 1956747 -252328.003
## Apr 2015 -844529.2119 1964016 -185249.556
## May 2015 -22366.3640 1966320 -73784.454
## Jun 2015 894329.5143 1968624 169915.618
## Jul 2015 1831172.5135 1970928 606454.569
## Aug 2015 2338719.8987 1959305 695438.764
## Sep 2015 1375035.0639 1947683 326982.179
## Oct 2015 11725.5715 1936060 -95106.749
## Nov 2015 -1036514.6333 1921706 -243732.880
## Dec 2015 -1134612.1812 1907351 -284838.668
## Jan 2016 -1126554.9520 1892996 -207381.233
## Feb 2016 -1206444.5114 1871952 -232319.173
## Mar 2016 -1102673.4276 1850907 -119512.756
## Apr 2016 -854581.2314 1829863 -73696.452
## May 2016 -23599.9988 1807609 168657.845
## Jun 2016 904106.7613 1785356 281533.615seasonal_component <- stl$time.series[,1]
ts_ds <- ts1 - seasonal_component
plot(ts_ds)
isSeasonal(ts_ds)## [1] FALSETest de la stationnarité:
Ce code utilise les fonctions ggAcf() et ggPacf() pour tracer la fonction d’autocorrélation (ACF) et la fonction d’autocorrélation partielle (PACF) de la série désaisonnalisée (ts_ds).
ggAcf(ts_ds) # ACF
ggPacf(ts_ds) # PACF
A première vue, on voit que celle-ci est clairement non stationnaire
puisque l’ACF décroît lentement au fil du temps, ce qui indique que les
observations sont fortement corrélées même à des distances temporelles
et importantes et la PACF montre une décroissance lente avec des retours
réguliers au-dessus du seuil de signification, ce qui indique que les
observations sont fortement corrélées même à des distances temporelles
importantes.
Nous pouvons en outre le confirmer avec le tests ADF : Le test Augmented Dickey-Fuller (ADF) est un test statistique utilisé pour tester la stationnarité d’une série temporelle. Il teste l’hypothèse nulle selon laquelle la série temporelle n’est pas stationnaire contre l’hypothèse alternative selon laquelle la série temporelle est stationnaire. Le test ADF peut être effectué de différentes manières, selon l’hypothèse alternative choisie. Dans ce cas, le test a été effectué avec trois hypothèses alternatives différentes : Type 1 : sans dérive ni tendance Type 2 : avec une dérive mais sans tendance Type 3 : avec une dérive et une tendance Le test ADF calcule une statistique de test appelée statistique ADF. Si cette statistique est inférieure à une certaine valeur critique, l’hypothèse nulle est rejetée, ce qui suggère que la série temporelle est stationnaire.
adf.test(ts_ds)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 0.0933 0.670
## [2,] 1 -0.3141 0.553
## [3,] 2 -0.1244 0.607
## [4,] 3 0.0805 0.666
## [5,] 4 0.1090 0.674
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -1.94 0.3531
## [2,] 1 -3.40 0.0143
## [3,] 2 -2.61 0.0963
## [4,] 3 -1.94 0.3520
## [5,] 4 -1.88 0.3771
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -3.07 0.1300
## [2,] 1 -5.16 0.0100
## [3,] 2 -4.32 0.0100
## [4,] 3 -3.55 0.0408
## [5,] 4 -3.59 0.0368
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Dans les résultats affichés, les valeurs de la statistique ADF et de la p-valeur pour chaque hypothèse alternative et chaque lag sont présentées. Si la p-valeur est inférieure au niveau de signification choisi (typiquement 0,05), l’hypothèse nulle est rejetée. Dans ce cas, les résultats montrent que pour les trois types d’hypothèses alternatives, la p-valeur est supérieure au niveau de signification de 0,05. Par conséquent, on ne peut pas rejeter l’hypothèse nulle et on peut conclure que la série temporelle est non stationnaire.
Stationnariser la série:
Ce code applique la fonction diff sur la série ts_ds avec une différence d’ordre 1, c’est-à-dire qu’il calcule la différence entre chaque observation et son observation précédente. Cela permet de transformer la série en une nouvelle série, où chaque observation représente la variation entre deux observations consécutives. Cette technique est souvent utilisée pour rendre une série temporelle stationnaire, car la différenciation peut aider à éliminer la tendance et à rendre la série plus aléatoire. La nouvelle série obtenue après différenciation est stockée dans l’objet ts_ds_diff.
ts_ds_diff<- diff(ts_ds, differences=1)
autoplot(ts_ds_diff) # plot
Nous pouvons exécuter le test ADF pour vérifier la stationnarité de la nouvelle série:
adf.test(ts_ds_diff)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -6.91 0.01
## [2,] 1 -7.45 0.01
## [3,] 2 -7.80 0.01
## [4,] 3 -6.37 0.01
## [5,] 4 -6.17 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -6.89 0.01
## [2,] 1 -7.43 0.01
## [3,] 2 -7.78 0.01
## [4,] 3 -6.37 0.01
## [5,] 4 -6.18 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -6.88 0.01
## [2,] 1 -7.42 0.01
## [3,] 2 -7.77 0.01
## [4,] 3 -6.36 0.01
## [5,] 4 -6.16 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Dans ce cas, les valeurs de ADF pour tous les trois types sont inférieures à la valeur critique de 1%, ce qui signifie que nous pouvons rejeter l’hypothèse nulle selon laquelle la série n’est pas stationnaire et accepter l’hypothèse alternative selon laquelle elle est stationnaire. Les p-values étant toutes inférieures ou égales à 0,01, cela confirme également que la série est effectivement stationnaire. De plus, on peut déduire qu’il y a absence de racines unitaires.
Modélisation avec ARIMA :
ARIMA (Autoregressive Integrated Moving Average) est un modèle classique de séries temporelles qui utilise l’autorégression, la moyenne mobile et la différenciation pour modéliser une série. Il peut être utilisé pour modéliser des séries stationnaires et non-stationnaires. D’abord, nous avons divisé notre série en deux, une partie train et une partie test
train1 <- window(ts_ds_diff, end=c(2014,12))
test1 <- window(ts_ds_diff,start=c(2015,1))La fonction auto.arima permet d’automatiquement sélectionner les ordres de l’ARIMA qui minimisent l’information de Bayes ou le critère d’information d’Akaike (AIC). Les arguments stepwise, approximation, seasonal et allowdrift sont utilisés pour spécifier les options de modélisation. stepwise=FALSE indique que la méthode stepwise de sélection de modèles ne sera pas utilisée et que l’algorithme de recherche complet sera exécuté. approximation=FALSE indique que la fonction d’approximation des erreurs ne sera pas utilisée. seasonal=FALSE indique que le modèle ne prendra pas en compte la saisonnalité. allowdrift=FALSE indique que la dérive n’est pas autorisée dans le modèle. L’argument trace=TRUE permet d’afficher les détails de l’ajustement de modèle en cours d’exécution.
ts_arima = auto.arima(train1,
stepwise=FALSE,
approximation=FALSE,
seasonal=FALSE,
allowdrift=FALSE,
trace=TRUE)##
## ARIMA(0,0,0) with zero mean : 2492.408
## ARIMA(0,0,0) with non-zero mean : 2494.421
## ARIMA(0,0,1) with zero mean : 2484.192
## ARIMA(0,0,1) with non-zero mean : 2486.282
## ARIMA(0,0,2) with zero mean : 2485.839
## ARIMA(0,0,2) with non-zero mean : 2487.981
## ARIMA(0,0,3) with zero mean : 2481.323
## ARIMA(0,0,3) with non-zero mean : 2482.657
## ARIMA(0,0,4) with zero mean : 2477.914
## ARIMA(0,0,4) with non-zero mean : 2478.949
## ARIMA(0,0,5) with zero mean : 2479.32
## ARIMA(0,0,5) with non-zero mean : 2480.283
## ARIMA(1,0,0) with zero mean : 2485.483
## ARIMA(1,0,0) with non-zero mean : 2487.58
## ARIMA(1,0,1) with zero mean : 2486.174
## ARIMA(1,0,1) with non-zero mean : 2488.314
## ARIMA(1,0,2) with zero mean : 2487.151
## ARIMA(1,0,2) with non-zero mean : 2489.338
## ARIMA(1,0,3) with zero mean : 2477.788
## ARIMA(1,0,3) with non-zero mean : 2478.709
## ARIMA(1,0,4) with zero mean : 2479.545
## ARIMA(1,0,4) with non-zero mean : 2480.535
## ARIMA(2,0,0) with zero mean : 2484.004
## ARIMA(2,0,0) with non-zero mean : 2486.119
## ARIMA(2,0,1) with zero mean : 2474.829
## ARIMA(2,0,1) with non-zero mean : 2475.65
## ARIMA(2,0,2) with zero mean : 2476.381
## ARIMA(2,0,2) with non-zero mean : 2477.18
## ARIMA(2,0,3) with zero mean : 2478.547
## ARIMA(2,0,3) with non-zero mean : 2479.458
## ARIMA(3,0,0) with zero mean : 2479.507
## ARIMA(3,0,0) with non-zero mean : 2481.501
## ARIMA(3,0,1) with zero mean : 2476.289
## ARIMA(3,0,1) with non-zero mean : 2477.105
## ARIMA(3,0,2) with zero mean : 2478.316
## ARIMA(3,0,2) with non-zero mean : 2479.201
## ARIMA(4,0,0) with zero mean : 2480.151
## ARIMA(4,0,0) with non-zero mean : 2482.015
## ARIMA(4,0,1) with zero mean : 2478.464
## ARIMA(4,0,1) with non-zero mean : 2479.357
## ARIMA(5,0,0) with zero mean : 2481.017
## ARIMA(5,0,0) with non-zero mean : 2482.684
##
##
##
## Best model: ARIMA(2,0,1) with zero meansummary(ts_arima)## Series: train1
## ARIMA(2,0,1) with zero mean
##
## Coefficients:
## ar1 ar2 ma1
## 1.0795 -0.4578 -0.8323
## s.e. 0.1072 0.0925 0.0829
##
## sigma^2 = 1.131e+10: log likelihood = -1233.19
## AIC=2474.38 AICc=2474.83 BIC=2484.6
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 11336.36 104645 70409.43 71.22572 205.9238 0.842597 -0.04335579La fonction auto.arima retourne un objet Arima qui contient les coefficients de l’ARIMA, les résidus, les statistiques du modèle et autres informations. Il est possible d’accéder à ces informations en utilisant les fonctions summary, coef, fitted, resid, etc. La sortie fournit les résultats de la meilleure modélisation de la série train1 obtenue à l’aide de la fonction auto.arima. Le modèle sélectionné est un ARIMA(2,0,1) avec une moyenne nulle, ce qui indique qu’il s’agit d’un modèle autorégressif intégré à la moyenne mobile d’ordre 2 et d’ordre de différenciation 0. Les coefficients estimés pour les termes AR et MA sont également fournis, ainsi que leurs erreurs standard correspondantes. Le carré de l’écart type de l’erreur de prédiction du modèle (sigma^2) est également donné, ainsi que le logarithme de la vraisemblance du modèle et les critères d’information tels que AIC, AICc et BIC. Les mesures d’erreur sur l’ensemble d’entraînement indiquent que le modèle a une erreur moyenne (ME) de 11336.36 et un écart type racine de la moyenne (RMSE) de 104645. Le MAE est la moyenne de la valeur absolue des erreurs de prédiction et est égal à 70409.43. Le MPE est l’erreur de prédiction en pourcentage moyen et le MAPE est l’erreur de prédiction absolue en pourcentage moyen. La MASE est une mesure d’erreur de prédiction relative par rapport à un modèle de marche aléatoire, tandis que la valeur de l’ACF1 est l’autocorrélation de la première différence de la série.
Validation du modèle ARIMA(2,0,1):
Il est recommandé de toujours inspecter visuellement les résidus du modèle pour s’assurer qu’ils sont bien des bruits blancs. On peut le faire en utilisant la fonction checkresiduals. Cette fonction produit plusieurs graphiques pour vérifier si les résidus sont indépendants, ont une distribution normale et si leur variance est constante.
checkresiduals(ts_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,1) with zero mean
## Q* = 10.672, df = 16, p-value = 0.8293
##
## Model df: 3. Total lags used: 19Ljung-Box test est utilisé pour tester l’autocorrélation des résidus d’un modèle. Le résultat montre que la valeur de Q* est 10.672 et le degré de liberté est 16. La p-valeur est 0.8293, ce qui indique qu’il n’y a pas suffisamment de preuves pour rejeter l’hypothèse nulle selon laquelle les résidus sont indépendants et identiquement distribués. Cela suggère que le modèle ARIMA(2,0,1) avec zéro moyen est adéquat pour modéliser la série temporelle.
Effectuons maintenant les prévisions: Le modèle ARIMA(2,0,1) est utilisé pour prédire les valeurs futures de la série à l’aide de la commande predict(), en spécifiant le nombre de pas de temps à prédire avec n.ahead = length(test1). Les prévisions sont stockées dans l’objet predictions. Le premier graphe affiche les prévisions générées par le modèle, tandis que le second graphe affiche la série de données de test test1. Enfin, la fonction accuracy() est utilisée pour calculer les mesures de précision de la prédiction en comparant les prévisions avec les valeurs réelles de la série de test test1. Les mesures de précision sont renvoyées à la sortie standard.
predictions <- predict(ts_arima, n.ahead = length(test1))$pred
ggplot() +
geom_line(aes(x = index(test1), y = coredata(test1), color = "Test")) +
geom_line(aes(x = index(predictions), y = predictions, color = "Forecast")) +
scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
ggtitle("Test et Forecast")
accuracy(predictions, test1)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 13430.73 201442.8 148301.4 -93.22257 309.4904 0.5686275 0.7097519Le modèle ARIMA semble donner des résultats relativement précis car les valeurs de RMSE, MAE et MAPE sont assez faibles. Cependant, le MPE est assez élevé, ce qui suggère que le modèle peut avoir une tendance à sous-estimer les valeurs. Il est également important de vérifier si le coefficient d’autocorrélation est significativement différent de zéro, car cela pourrait indiquer des erreurs systématiques dans les prévisions. Le Theil’s U est inférieur à 1, ce qui indique que le modèle est meilleur que la prévision naïve.
Modélisation avec SARIMA:
SARIMA (Seasonal ARIMA) est une extension de ARIMA pour les séries saisonnières. Il utilise l’autorégression, la moyenne mobile, la différenciation saisonnière et la différenciation normale pour modéliser une série saisonnière. Commençons d’abord par stationnariser la série saisonnière ts_ds:
ts1_stat <- diff(ts1, 1)
adf.test(ts1_stat)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -4.40 0.01
## [2,] 1 -8.08 0.01
## [3,] 2 -7.05 0.01
## [4,] 3 -5.69 0.01
## [5,] 4 -6.42 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -4.39 0.01
## [2,] 1 -8.04 0.01
## [3,] 2 -7.02 0.01
## [4,] 3 -5.66 0.01
## [5,] 4 -6.39 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -4.36 0.01
## [2,] 1 -8.00 0.01
## [3,] 2 -6.98 0.01
## [4,] 3 -5.62 0.01
## [5,] 4 -6.35 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Les résultats présentés ici montrent trois types de tests ADF avec différents niveaux de tendance ou de dérive. Dans les trois types, les valeurs de la statistique ADF (colonne ADF) sont inférieures aux valeurs critiques à un niveau de signification de 0,05, ce qui suggère que la série temporelle est stationnaire. De plus, les valeurs de p (colonne p.value) sont toutes inférieures à 0,05, ce qui indique que la probabilité d’observer ces résultats par hasard est très faible.
Vérifions la saisonnalité de la nouvelle série: ce test a été défini précédemment.
isSeasonal(ts1_stat)## [1] TRUESélectionner la période d’entraînement de la série temporelle ts1_stat allant de janvier 2006 à décembre 2014 et stocker cette période dans l’objet train2. Sélectionner la période de test de la série temporelle ts1_stat allant de janvier 2015 jusqu’à la fin de la série et stocker cette période dans l’objet test2.
train2 <- window(ts1_stat, end=c(2014,12))
test2<- window(ts1_stat,start=c(2015,1))Appliquons la fonction auto.arima() sur la série temporelle train2 pour trouver le meilleur modèle SARIMA. Les arguments stepwise=FALSE et approximation=FALSE sont utilisés pour que la fonction teste tous les modèles possibles au lieu de se limiter aux modèles les plus simples. L’argument seasonal=TRUE est utilisé pour spécifier que la série est saisonnière, et l’argument allowdrift=FALSE est utilisé pour exclure les modèles avec une dérive. L’argument trace=TRUE est utilisé pour afficher les informations de progression pendant l’exécution de la fonction. Puis, on affiche le résumé du modèle SARIMA sélectionné avec la fonction summary().
ts_sarima = auto.arima(train2,
stepwise=FALSE,
approximation=FALSE,
seasonal=TRUE,
allowdrift=FALSE,
trace=TRUE)##
## ARIMA(0,0,0)(0,1,0)[12] : 2196.707
## ARIMA(0,0,0)(0,1,1)[12] : 2198.693
## ARIMA(0,0,0)(0,1,2)[12] : 2200.839
## ARIMA(0,0,0)(1,1,0)[12] : 2198.695
## ARIMA(0,0,0)(1,1,1)[12] : Inf
## ARIMA(0,0,0)(1,1,2)[12] : 2203.024
## ARIMA(0,0,0)(2,1,0)[12] : 2200.846
## ARIMA(0,0,0)(2,1,1)[12] : Inf
## ARIMA(0,0,0)(2,1,2)[12] : Inf
## ARIMA(0,0,1)(0,1,0)[12] : 2197.837
## ARIMA(0,0,1)(0,1,1)[12] : 2199.964
## ARIMA(0,0,1)(0,1,2)[12] : 2202.104
## ARIMA(0,0,1)(1,1,0)[12] : 2199.967
## ARIMA(0,0,1)(1,1,1)[12] : Inf
## ARIMA(0,0,1)(1,1,2)[12] : Inf
## ARIMA(0,0,1)(2,1,0)[12] : 2202.05
## ARIMA(0,0,1)(2,1,1)[12] : 2204.278
## ARIMA(0,0,1)(2,1,2)[12] : Inf
## ARIMA(0,0,2)(0,1,0)[12] : 2199.304
## ARIMA(0,0,2)(0,1,1)[12] : 2201.491
## ARIMA(0,0,2)(0,1,2)[12] : 2203.623
## ARIMA(0,0,2)(1,1,0)[12] : 2201.494
## ARIMA(0,0,2)(1,1,1)[12] : Inf
## ARIMA(0,0,2)(1,1,2)[12] : Inf
## ARIMA(0,0,2)(2,1,0)[12] : 2203.597
## ARIMA(0,0,2)(2,1,1)[12] : Inf
## ARIMA(0,0,3)(0,1,0)[12] : 2197.078
## ARIMA(0,0,3)(0,1,1)[12] : 2199.304
## ARIMA(0,0,3)(0,1,2)[12] : 2201.296
## ARIMA(0,0,3)(1,1,0)[12] : 2199.312
## ARIMA(0,0,3)(1,1,1)[12] : Inf
## ARIMA(0,0,3)(2,1,0)[12] : 2201.09
## ARIMA(0,0,4)(0,1,0)[12] : 2196.239
## ARIMA(0,0,4)(0,1,1)[12] : 2198.316
## ARIMA(0,0,4)(1,1,0)[12] : 2198.36
## ARIMA(0,0,5)(0,1,0)[12] : 2196.748
## ARIMA(1,0,0)(0,1,0)[12] : 2198.02
## ARIMA(1,0,0)(0,1,1)[12] : 2200.167
## ARIMA(1,0,0)(0,1,2)[12] : 2202.335
## ARIMA(1,0,0)(1,1,0)[12] : 2200.168
## ARIMA(1,0,0)(1,1,1)[12] : Inf
## ARIMA(1,0,0)(1,1,2)[12] : Inf
## ARIMA(1,0,0)(2,1,0)[12] : 2202.309
## ARIMA(1,0,0)(2,1,1)[12] : 2204.538
## ARIMA(1,0,0)(2,1,2)[12] : Inf
## ARIMA(1,0,1)(0,1,0)[12] : 2199.769
## ARIMA(1,0,1)(0,1,1)[12] : 2201.972
## ARIMA(1,0,1)(0,1,2)[12] : 2204.146
## ARIMA(1,0,1)(1,1,0)[12] : 2201.972
## ARIMA(1,0,1)(1,1,1)[12] : Inf
## ARIMA(1,0,1)(1,1,2)[12] : Inf
## ARIMA(1,0,1)(2,1,0)[12] : 2204.09
## ARIMA(1,0,1)(2,1,1)[12] : 2206.357
## ARIMA(1,0,2)(0,1,0)[12] : 2193.337
## ARIMA(1,0,2)(0,1,1)[12] : 2195.599
## ARIMA(1,0,2)(0,1,2)[12] : 2197.705
## ARIMA(1,0,2)(1,1,0)[12] : 2195.6
## ARIMA(1,0,2)(1,1,1)[12] : Inf
## ARIMA(1,0,2)(2,1,0)[12] : 2197.591
## ARIMA(1,0,3)(0,1,0)[12] : 2195.099
## ARIMA(1,0,3)(0,1,1)[12] : 2197.301
## ARIMA(1,0,3)(1,1,0)[12] : 2197.322
## ARIMA(1,0,4)(0,1,0)[12] : 2197.144
## ARIMA(2,0,0)(0,1,0)[12] : 2199.146
## ARIMA(2,0,0)(0,1,1)[12] : 2201.353
## ARIMA(2,0,0)(0,1,2)[12] : 2203.46
## ARIMA(2,0,0)(1,1,0)[12] : 2201.354
## ARIMA(2,0,0)(1,1,1)[12] : Inf
## ARIMA(2,0,0)(1,1,2)[12] : Inf
## ARIMA(2,0,0)(2,1,0)[12] : 2203.375
## ARIMA(2,0,0)(2,1,1)[12] : 2205.586
## ARIMA(2,0,1)(0,1,0)[12] : 2192.584
## ARIMA(2,0,1)(0,1,1)[12] : 2194.643
## ARIMA(2,0,1)(0,1,2)[12] : 2196.703
## ARIMA(2,0,1)(1,1,0)[12] : 2194.687
## ARIMA(2,0,1)(1,1,1)[12] : 2196.904
## ARIMA(2,0,1)(2,1,0)[12] : 2196.438
## ARIMA(2,0,2)(0,1,0)[12] : 2193.279
## ARIMA(2,0,2)(0,1,1)[12] : 2194.27
## ARIMA(2,0,2)(1,1,0)[12] : 2194.498
## ARIMA(2,0,3)(0,1,0)[12] : Inf
## ARIMA(3,0,0)(0,1,0)[12] : 2200.369
## ARIMA(3,0,0)(0,1,1)[12] : 2202.501
## ARIMA(3,0,0)(0,1,2)[12] : 2204.704
## ARIMA(3,0,0)(1,1,0)[12] : 2202.523
## ARIMA(3,0,0)(1,1,1)[12] : 2204.797
## ARIMA(3,0,0)(2,1,0)[12] : 2204.583
## ARIMA(3,0,1)(0,1,0)[12] : 2194.479
## ARIMA(3,0,1)(0,1,1)[12] : 2196.395
## ARIMA(3,0,1)(1,1,0)[12] : 2196.475
## ARIMA(3,0,2)(0,1,0)[12] : 2195.516
## ARIMA(4,0,0)(0,1,0)[12] : 2200.986
## ARIMA(4,0,0)(0,1,1)[12] : 2202.995
## ARIMA(4,0,0)(1,1,0)[12] : 2203.037
## ARIMA(4,0,1)(0,1,0)[12] : 2196.276
## ARIMA(5,0,0)(0,1,0)[12] : 2200.956
##
##
##
## Best model: ARIMA(2,0,1)(0,1,0)[12]summary(ts_sarima)## Series: train2
## ARIMA(2,0,1)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ma1
## 0.9312 -0.2398 -0.9309
## s.e. 0.1165 0.1100 0.0606
##
## sigma^2 = 1.612e+10: log likelihood = -1092.04
## AIC=2192.07 AICc=2192.58 BIC=2201.75
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 11056.02 116507.2 71381.5 -30.05839 78.53752 0.8132461 -0.02493748Le meilleur modèle sélectionné est ARIMA(2,0,1)(0,1,0)[12]. Cela signifie qu’il s’agit d’un modèle autorégressif intégré de moyenne mobile saisonnier d’ordre 2, 0 et 1 pour les termes AR, I et MA, respectivement. Le modèle inclut également une différence saisonnière d’ordre 1 avec une période de 12 (saison). Les coefficients associés aux termes AR1, AR2 et MA1 sont respectivement de 0,9312, -0,2398 et -0,9309. La valeur p associée à chaque coefficient est inférieure à 0,05, ce qui suggère que les coefficients sont significatifs. La valeur de sigma² est de 1,612e+10 et la log-vraisemblance associée au modèle est de -1092,04. Les mesures d’erreur d’entraînement (ME, RMSE, MAE, MPE, MAPE, MASE et ACF1) sont également fournies, ce qui permet d’évaluer la performance du modèle sur l’ensemble d’entraînement.
Validation du modèle ARIMA(2,0,1)(0,1,0)[12]:
Validation par résidus : cette méthode consiste à vérifier si les résidus du modèle ARIMA sont aléatoires et ne présentent pas de structure temporelle. Pour cela, on peut utiliser les fonctions checkresiduals ou tsdiag, qui tracent différents graphiques pour vérifier la normalité, l’absence d’autocorrélation, l’homoscédasticité, etc.
checkresiduals(ts_sarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,1)(0,1,0)[12]
## Q* = 15.529, df = 16, p-value = 0.4863
##
## Model df: 3. Total lags used: 19Le test de Ljung-Box est utilisé pour vérifier si les résidus d’un modèle ARIMA sont aléatoires ou s’ils présentent une autocorrélation significative. Dans ce cas, le modèle ARIMA sélectionné est (2,0,1)(0,1,0)[12]. Le résultat du test de Ljung-Box montre un Q* de 15.529 et un degré de liberté de 16, avec une p-valeur de 0.4863. Comme la p-valeur est supérieure à 0,05 (niveau de signification communément choisi), on ne peut pas rejeter l’hypothèse nulle selon laquelle les résidus sont aléatoires. Cela indique que le modèle ARIMA(2,0,1)(0,1,0)[12] est un bon ajustement pour les données.
Effectuons maintenant les prévisions:
predictions2 <- predict(ts_sarima, n.ahead = length(test2))$pred
ggplot() +
geom_line(aes(x = index(test2), y = coredata(test2), color = "Test")) +
geom_line(aes(x = index(predictions2), y = predictions2, color = "Forecast")) +
scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
ggtitle("Test et Forecast")
accuracy(predictions2, test2)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -11429.78 72154.85 58477.91 15.59264 24.66167 -0.2513726 0.2228339Le tableau affiche les erreurs de prévision du modèle SARIMA pour la période de test. Les différentes métriques d’erreur sont les suivantes : Les valeurs de RMSE et MAE sont relativement élevées, ce qui indique que les erreurs de prévision ne sont pas négligeables. La valeur de MPE est positive, ce qui signifie que les prévisions sont en moyenne supérieures aux valeurs réelles. La valeur de MAPE est inférieure à 25%, ce qui indique que les prévisions ne sont pas très éloignées des valeurs réelles. La valeur de l’ACF1 est négative, ce qui suggère que les erreurs de prévision ont tendance à s’annuler mutuellement. La valeur de Theil’s U est relativement faible, ce qui indique que le modèle n’est pas très précis dans ses prévisions.
Modélisation par la méthode de Regression:
La méthode de régression est une méthode couramment utilisée pour la modélisation des séries temporelles. Elle consiste à modéliser la série temporelle en tant que fonction d’autres variables explicatives, telles que le temps, des indicateurs économiques ou météorologiques, ou d’autres facteurs pertinents. Pour appliquer la modélisation par régression linéaire sur une série temporelle, il est généralement préférable que la série temporelle soit stationnaire. La stationnarité d’une série temporelle signifie que la moyenne, la variance et la covariance ne dépendent pas du temps et restent constantes sur la période étudiée. Pour celà, nous allons appliquer la regression sur la série “ts1_stat”.
On a dévisé la série en deux sous séries train et test:
train3 <- window(ts1_stat, end=c(2014,12))
test3<- window(ts1_stat,start=c(2015,1))La ligne de code ts_reg <- tslm(ts1_stat ~ trend + season(ts1_stat)) effectue la modélisation d’une série temporelle ts1_stat à l’aide d’un modèle de régression linéaire multiple avec une tendance linéaire et une composante saisonnière. Plus précisément, la variable dépendante est ts1_stat et les variables indépendantes sont trend et season(ts1_stat). La fonction season() est utilisée pour extraire la composante saisonnière de la série temporelle. La fonction tslm() est une fonction qui permet d’estimer les paramètres d’un modèle de régression linéaire multiple. La ligne de code suivante summary(ts_reg) affiche un résumé des résultats de la modélisation, qui incluent des informations sur les coefficients de régression, les statistiques de diagnostic, les erreurs de prévision et la qualité globale du modèle.
ts_reg <- tslm(ts1_stat ~ trend + season(ts1_stat))
summary(ts_reg)##
## Call:
## tslm(formula = ts1_stat ~ trend + season(ts1_stat))
##
## Residuals:
## Min 1Q Median 3Q Max
## -568717 -61364 8645 68317 596552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -28082.5 67514.0 -0.416 0.67834
## trend 307.5 519.0 0.592 0.55492
## season(ts1_stat)February -39656.5 82610.8 -0.480 0.63225
## season(ts1_stat)March 117982.1 82596.1 1.428 0.15629
## season(ts1_stat)April 261784.5 82584.7 3.170 0.00202 **
## season(ts1_stat)May 795354.9 82576.5 9.632 6.36e-16 ***
## season(ts1_stat)June 829335.8 82571.7 10.044 < 2e-16 ***
## season(ts1_stat)July 929261.6 84772.3 10.962 < 2e-16 ***
## season(ts1_stat)August 513004.5 84754.8 6.053 2.50e-08 ***
## season(ts1_stat)September -896899.4 84740.5 -10.584 < 2e-16 ***
## season(ts1_stat)October -1234891.3 84729.4 -14.575 < 2e-16 ***
## season(ts1_stat)November -996793.7 84721.4 -11.766 < 2e-16 ***
## season(ts1_stat)December -73445.3 84716.6 -0.867 0.38804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 179700 on 100 degrees of freedom
## Multiple R-squared: 0.9436, Adjusted R-squared: 0.9368
## F-statistic: 139.4 on 12 and 100 DF, p-value: < 2.2e-16Le résumé du modèle donne une liste des coefficients estimés pour chaque variable explicative incluse dans le modèle, ainsi que leurs erreurs standards, t-values et p-values. Le coefficient de détermination (R2) multiple est également fourni, qui est une mesure de la qualité de l’ajustement global du modèle. Dans ce cas, le R2 multiple est égal à 0,9436, ce qui indique que le modèle ajusté explique environ 94,36% de la variation dans les données. Le modèle indique que toutes les variables saisonnières (mois) sont significatives avec des p-values inférieures à 0,05, ce qui suggère qu’elles ont une influence significative sur la série temporelle. Cependant, le coefficient d’interception et le coefficient de tendance ne sont pas significatifs avec des p-values supérieures à 0,05. Cela signifie qu’il n’y a pas suffisamment de preuves pour conclure que la série temporelle a une tendance significative ou un niveau significatif différent de zéro.
Affichons maintenant le vecteur contenant les valeurs prédites par le modèle de régression pour chaque observation de la série temporelle ts1_stat. Il s’agit des valeurs qui minimisent l’erreur de prédiction entre les observations réelles et les prédictions du modèle.
ts_reg$fitted.values## Jan Feb Mar Apr May
## 2007 -67431.570 90514.530 234624.430 768502.230
## 2008 -24392.862 -63741.910 94204.190 238314.090 772191.890
## 2009 -20703.202 -60052.250 97893.850 242003.750 775881.550
## 2010 -17013.542 -56362.590 101583.510 245693.410 779571.210
## 2011 -13323.882 -52672.930 105273.170 249383.070 783260.870
## 2012 -9634.222 -48983.270 108962.830 253072.730 786950.530
## 2013 -5944.562 -45293.610 112652.490 256762.390 790640.190
## 2014 -2254.902 -41603.950 116342.150 260452.050 794329.850
## 2015 1434.758 -37914.290 120031.810 264141.710 798019.510
## 2016 5124.418 -34224.630 123721.470 267831.370 801709.170
## Jun Jul Aug Sep Oct
## 2007 802790.630 903023.915 487074.249 -922522.196 -1260206.529
## 2008 806480.290 906713.575 490763.909 -918832.536 -1256516.869
## 2009 810169.950 910403.235 494453.569 -915142.876 -1252827.209
## 2010 813859.610 914092.895 498143.229 -911453.216 -1249137.549
## 2011 817549.270 917782.556 501832.889 -907763.556 -1245447.889
## 2012 821238.930 921472.216 505522.549 -904073.895 -1241758.229
## 2013 824928.590 925161.876 509212.209 -900384.235 -1238068.569
## 2014 828618.250 928851.536 512901.869 -896694.575 -1234378.909
## 2015 832307.910 932541.196 516591.529 -893004.915 -1230689.249
## 2016 835997.570
## Nov Dec
## 2007 -1021801.529 -98145.640
## 2008 -1018111.869 -94455.980
## 2009 -1014422.209 -90766.320
## 2010 -1010732.549 -87076.660
## 2011 -1007042.889 -83387.000
## 2012 -1003353.229 -79697.340
## 2013 -999663.569 -76007.680
## 2014 -995973.909 -72318.020
## 2015 -992284.249 -68628.360
## 2016autoplot(ts_reg$fitted.values)
Validation de la méthode de regression:
checkresiduals(ts_reg)
##
## Breusch-Godfrey test for serial correlation of order up to 23
##
## data: Residuals from Linear regression model
## LM test = 78.478, df = 23, p-value = 5.59e-08Le test de Breusch-Godfrey est utilisé pour détecter la présence de corrélation série temporelle dans les résidus d’un modèle de régression. Dans ce cas, le test a été effectué sur les résidus d’un modèle de régression linéaire. La statistique de test LM est de 78,478 avec 23 degrés de liberté, et la valeur de p est de 5,59e-08. Puisque la valeur de p est inférieure à 0,05, on rejette l’hypothèse nulle selon laquelle il n’y a pas de corrélation série temporelle dans les résidus, et on conclut qu’il y a une corrélation série temporelle dans les résidus. Cela suggère que le modèle de régression linéaire n’est pas suffisant pour modéliser la série temporelle.
Modélisation par la méthode de lissage exponentielle:
Les modèles de lissage exponentiel sont des méthodes de prévision de séries temporelles qui s’appuient sur la notion de moyenne mobile exponentielle pour estimer les tendances locales et saisonnières d’une série temporelle. Ces modèles ne sont pas linéaires, car ils utilisent des coefficients de lissage exponentiels pour donner plus ou moins d’importance aux observations passées en fonction de leur éloignement dans le temps. Les modèles de lissage exponentiel peuvent être appliqués dans des situations où la série temporelle a une tendance, une saisonnalité et des composantes d’erreur qui sont stables ou qui changent lentement au fil du temps. En d’autres termes, ils sont adaptés aux séries temporelles qui ne présentent pas de tendances ou de motifs de saisonnalité complexes et qui ne présentent pas de changements brusques dans les modèles de comportement. Le lissage exponentiel est souvent utilisé pour modéliser des séries temporelles non stationnaires et saisonnières. La série doit être lissée pour éliminer les variations de courte durée et mettre en évidence les tendances sous-jacentes qui peuvent être utilisées pour la prévision. Les modèles de lissage exponentiel sont donc appliqués directement sur la série temporelle originale.
On divise d’abord la série originale en deux sous séries:
train4 <- window(ts1, end=c(2014,12))
test4 <- window(ts1,start=c(2015,1))La commande ts1_hw <- HoltWinters(train4, beta=TRUE, gamma=TRUE) applique le modèle de Holt-Winters avec tendance et saisonnalité additive sur la série train4 pour estimer les paramètres nécessaires à la prévision des valeurs futures. La tendance est modélisée par l’exponentielle double, et la saisonnalité est prise en compte par la méthode de lissage exponentiel triple. Les paramètres beta et gamma permettent de déterminer si la tendance et/ou la saisonnalité sont modélisées de manière additive ou multiplicative.
ts1_hw <- HoltWinters(train4, beta=TRUE, gamma=TRUE)
predictions4 <- predict(ts1_hw, n.ahead=18)
ggplot() +
geom_line(aes(x = index(test4), y = coredata(test4), color = "Test")) +
geom_line(aes(x = index(predictions4), y = predictions4, color = "Forecast")) +
scale_color_manual(values = c("Test" = "blue", "Forecast" = "darkred")) +
labs(x = "t", y = "Yt:Nombre de touristes", color = "Legend") +
ggtitle("Test et Forecast")
accuracy(predictions4, test4)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 1005335 1170533 1005335 78.99331 78.99331 0.7117418 2.075206ME : Erreur moyenne, c’est la moyenne des différences entre les valeurs observées et prédites. Ici, elle vaut 1005335, ce qui signifie que les prédictions ont tendance à sous-estimer les valeurs observées. RMSE : Erreur quadratique moyenne, c’est la racine carrée de la moyenne des carrés des différences entre les valeurs observées et prédites. Ici, elle vaut 1170533, ce qui représente la déviation standard des erreurs entre les valeurs observées et prédites. MAE : Erreur absolue moyenne, c’est la moyenne des valeurs absolues des différences entre les valeurs observées et prédites. Ici, elle vaut 1005335, ce qui indique une tendance de sous-estimation des prédictions. MPE : Erreur de pourcentage moyen, c’est la moyenne des différences en pourcentage entre les valeurs observées et prédites. Ici, elle vaut 78.99331, ce qui indique que les prédictions sont en moyenne inférieures de 78,99% aux valeurs observées. MAPE : Erreur absolue de pourcentage moyen, c’est la moyenne des valeurs absolues des différences en pourcentage entre les valeurs observées et prédites. Ici, elle vaut également 78.99331, ce qui confirme la tendance de sous-estimation des prédictions. ACF1 : Le coefficient de corrélation autocorrélatif de la première différence. S’il est proche de 1 ou -1, cela indique une forte autocorrélation dans les données. Ici, il vaut 0.7117418, ce qui indique une certaine autocorrélation dans les résidus du modèle. Theil’s U : Coefficient d’efficacité de la prédiction, il mesure l’efficacité de la prédiction par rapport à une méthode de référence telle que la prédiction naïve. Ici, il vaut 2.075206, ce qui indique que le modèle de lissage exponentiel est deux fois plus efficace que la méthode de prédiction naïve.
Comparaison des 3 modèles:
ARIMA_res = data.frame(accuracy(predictions, test1))
SARIMA_res = data.frame(accuracy(predictions2, test2))
LISSAGE_res = data.frame(accuracy(predictions4, test4))
resf <- data.frame(MAE=c(ARIMA_res$MAE,SARIMA_res$MAE,LISSAGE_res$MAE),MAPE = c(ARIMA_res$MAPE, SARIMA_res$MAPE, LISSAGE_res$MAPE),RMSE=c(ARIMA_res$RMSE,SARIMA_res$RMSE,LISSAGE_res$RMSE))
rownames(resf) <- c("ARIMA", "SARIMA", "LISSAGE EXPONENTIELLE")
resf## MAE MAPE RMSE
## ARIMA 148301.44 309.49035 201442.75
## SARIMA 58477.91 24.66167 72154.85
## LISSAGE EXPONENTIELLE 1005335.43 78.99331 1170533.25==> Dans ce cas, le modèle SARIMA a les valeurs les plus faibles pour les trois mesures, ce qui indique qu’il a la meilleure performance de prévision parmi les trois modèles testés. Le modèle de lissage exponentiel a les valeurs les plus élevées pour les trois mesures, ce qui indique qu’il a la pire performance de prévision.
PDF <- predict(ts_sarima, n.ahead=48)
orig_series <- cumsum(PDF$pred)
orig_series1 <- ts(orig_series,start=c(2015,1), end=c(2018,12), frequency = 12)
autoplot(orig_series1)